C.................... CALODDN.FOR;8 ..........6-JUN-1994 10:33:04.41 
C---- This subroutine is the main sequencing control for the odd N
C---- and vibrationally excited N2 solutions on the vertical grid. 
C---- It was developed in 1993 by splitting up the old RSVLPE.FOR
C---- into CALODDN and CALVN2. It is responsible for calling the
C---- odd N and N2* solutions in the proper order.
C---- Written by Phil Richards in February 1993.
C.... CALDIF, OFIJ, and VFIJ modified by P. Richards in April 2003 
C.... to include eddy diffusion.

      SUBROUTINE VERCON(J250,NOCON,S,JMAX,IDIM,DTV,VBLAT,GRD,ISOL,IVLPS
     >   ,IVIBN2,HMF2N,HMF2S,SATN,SATF,CHIN,CHIF)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      INTEGER NOCON,J250,J250C
      REAL SATN,SATF,CHIN,CHIF,GRD
      DIMENSION S(IDIM,1)       ! f90 added

      NOCON=0      !.... flag for non-convergence
C----- Northern hemisphere odd nitrogen
      CALL CALODD(J250,NOCON,1,JMAX/2,S,IDIM,DTV,VBLAT,GRD,-1,ISOL,17
     >   ,IVLPS)

C----- Northern hemisphere Vibrational N2 DENSITIES
      IF(IVIBN2.GE.1) 
     >  CALL CALVN2(1,JMAX/2,S,999,DTV,VBLAT,GRD,+1,ISOL,17,IVLPS
     >     ,HMF2N,HMF2S,SATN,SATF,CHIN,CHIF)

      J250C=JMAX+1-J250
C----- Southern hemisphere odd nitrogen
      CALL CALODD(J250C,NOCON,(JMAX+2)/2,JMAX,S,IDIM,DTV,-VBLAT,GRD,-1
     >   ,ISOL,17,IVLPS)

C----- Southern hemisphere Vibrational N2 DENSITIES
      IF(IVIBN2.GE.1) 
     >  CALL CALVN2((JMAX+2)/2,JMAX,S,999,DTV,-VBLAT,GRD,+1
     >      ,ISOL,17,IVLPS,HMF2N,HMF2S,SATN,SATF,CHIN,CHIF)

      RETURN
      END
C:::::::::::::::::::::::::::::::: CALODD ::::::::::::::::::::::::
C------ This subroutine is the main sequencing control for the odd N
C------ solutions. It was developed in 1993 by splitting up the old
C------ RSVLPE.FOR. It specifies time and convergence looping and calls
C------ subroutines to set up and solve the coupled continuity and 
C------ momentum equations using Newton's method.
C------ Written by Phil Richards in February 1993.
C------ Modified by P. Richards in December 2008 to allow for variable
C------ upper boundary altitude that is needed for low L-values. 
C------ Vertical grid cannot exceed the L-shell apex altitude so that 
C------ ion densities can be interpolated to vertical grid
      SUBROUTINE CALODD(J250,NOCON,MINJ,MAXJ,S,IDIM,DTV,VBLAT,GRD,IC,ISOL
     >   ,MAXIT,IVLPS)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      INTEGER NION,NEQ,NFLAG,NOCON,VD,LAXIT,J250
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      REAL AP,DEC,ETRAN,BLON,F107,F107A,SEC,GRD
      DIMENSION S(IDIM,1),PROD(6),ALOSS(6),VCON(VD),AM(7)
      COMMON/SLV/DELTA(1806),RHS(1806),WORK(50000)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/FJS/N(4,401),TI(3,401),FF(20)
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      COMMON/ALT/Z(401),DT,DW,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,NION
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/LPS/SZA(401),D1,D2,ISPC,I1,I2,IPRIN,JTIMAX,IPP,ISKP
      COMMON/VAGRID/JBV,JVB,JNB
C
       DATA  EPSN, RE  , DCR  ,DC, EPCON , NFLAG, IPR
     >   /    .8 ,6357., 0.0  ,.5, 1.0E-2,  0   ,  3  /
C------ Atomic masses for the odd N species, vibrational N2 factors
      DATA AM/14.,14.,30.,16.,3*14./,         VCON/VD*0.0D0/
C------- ID= unit number for output file.
        ID=3
        IF(MINJ.GT.1) ID=4
C------- setting up # grid points and time step
        JVMAX=JMAX
        VDT=DTV
        JTER=0

      !.. Find the upper boundary for vertical drift, making sure it is 
      !.. not greater than the field line apex altitude
      ! Modified in December 2008 to allow for variable upper boundary 
      DO J=1,JBV
        IF(MAXJ.LT.(JMAX+1)/2.AND.ZR(J).LT.Z(MAXJ)) JBP1=J
        IF(MAXJ.GT.(JMAX+1)/2.AND.ZR(J).LT.Z(MINJ)) JBP1=J
      ENDDO
      JBP=JBP1-1

      JBL=JVB           !.. Lower boundary index
      MIT=JBP-JBL+1     !.. # points on the altiude grid to solve
      NEQ=IPR*(MIT-1)   !.. # of equations to solve

C------- set up the neutral atmosphere on vertical grid
      CALL VATMOS(J250,MINJ,MAXJ,GRD,VBLAT,IC,ISOL,MAXIT,IVLPS,JBP)
C
C-------- Set up terms in the DE that do not change during iteration
      CALL CALDIF(JBP1,IPR,AM,IC)
C
C------ Fill the solution arrays NV with initial values or read previous 
C------ values from STR(IS,J). JCON stipulates northern or southern 
C------ hemisphere for storing in STR. VCON is N2* factor
      DO 21 J=1,JBP1
         NV(1,J)=N2D(J)
         NV(2,J)=N4S(J)
         NV(3,J)=NNO(J)
         NV(4,J)=O1D(J)
         JCON=J
         IF(MINJ.GT.1) JCON=MAXJ+1-J
         IF(STR(7,JCON).LT.1.0.OR.STR(7,JCON).GT.5) STR(7,JCON)=1.0
         VCON(J)=STR(7,JCON)
         IF(ISOL.EQ.0) CALL PRLN(1,0,J,0,IPR,PROD,ALOSS,VCON,IVLPS)
         DO 25 I=1,IPR
            IF(ISOL.NE.0) NV(I,J)=STR(I+7,JCON)
            IF(ISOL.NE.0) GO TO 25
C....            IF(I.NE.2) NV(I,J)=PROD(I)/ALOSS(I)
 25    CONTINUE
 21    CONTINUE
C
C------- save density at previous time step for calculation of dn/dt
      DO 41 J=1,JBP1
      DO 41 IS=1,IPR
         NVS(IS,J)=NV(IS,J)
 41    CONTINUE
C
C-------- set lower boundary conditions
      DO 99 J=1,JBL
C------- safety feature to prevent runaway
         IF(NV(2,J).GT.NV(3,J)) THEN
            NV(2,J)=1.0E4
            NV(3,J)=3.0E7
         ENDIF
C------- Iterate on previous values to improve N(4S) and NO
         DO 69 ITYMS=1,2
            CALL PRLN(1,0,J,0,IPR,PROD,ALOSS,VCON,IVLPS)
            DO 68 IS=1,IPR
              NV(IS,J)=PROD(IS)/ALOSS(IS)
 68         CONTINUE
 69      CONTINUE
 99   CONTINUE
C
      IF(ISOL.EQ.9) RETURN
C
C------- MAIN loop - begin iterative solution procedure ends at next row of
C**********************************************************************
      LAXIT=IABS(MAXIT)
      DO 220 JTER=1,LAXIT
C
C----- compute Fij profile and S matrix which BDSLV uses to
C----- solve the linear system ((S))*(DELTA)= (Fij)
      DO 145 J=2,MIT
         JC=J+JBL-1
         CALL OFIJ(JC,0,IPR,JBP,IC)
         KR=IPR*(J-2)
         DO 139 IRHS=1,IPR
 139     RHS(KR+IRHS)=F(IRHS)
 145  CONTINUE
C----- IJM is a switch for calling OMATRIX or not ...Cuts down on the amount
C----- of work in building the matrix.
       IJM=0
       IF((JTER/3)*3.EQ.JTER) IJM=1
       IF(JTER.LE.3) IJM=1
       IF(DCR.NE.1) IJM=1
C
       IF(IJM.EQ.1) CALL OMATRX(S,RHS,NEQ,IPR,MIT,JBL,JBP,IC)
C
C------- invert the Jacobian matrix S and record error comments
C-------   from the inversion routine BDSLV. IBW = band width
       IBW=2*IPR-1
       CALL BDSLV(NEQ,IBW,S,0,RHS,DELTA,WORK,NFLAG)
       IF(NFLAG.NE.0) THEN
         WRITE(6,*) ' ** ERROR Convergence problem in BDSLV *******'
         WRITE(6,*) ' Problem in odd N solution'
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP 3
      ENDIF
C
C----- Update previous densities NV with increment DELTA from BDSLV
C----- Test for convergence when IDIV=0 ... no densities change by more
C----- than specified amount.
 55    IDIV=0
       DCR=1.0
C------- First, test delta to see if it is bigger than NV -> -ve densities
C------- Set DCRP to ensure NV remains positive
      DO 142 I=1,IPR
         ION=2*IPR-I
         DO 142 J=2,MIT
         DCRP=1.0
         DINC=DELTA(IPR*J-ION)
         JC=J+JBL-1
         IF(DINC.LT.0.0) GO TO 142
         IF(ABS(DINC/NV(I,JC)).GT.EPSN) DCRP=DC*ABS(NV(I,JC)/DINC)
         IF(DCRP.LT.DCR) DCR=DCRP
 142  CONTINUE
C
C----- subtract all or part of DELTA from density NV to get new density
C----- for next iteration
 144  DO 42 J=2,MIT
         JC=J+JBL-1
         DO 42 I=1,IPR
         ION=2*IPR-I
         DINC=DELTA(IPR*J-ION)
         NV(I,JC)=NV(I,JC)-DINC*DCR
         IF(ABS(DINC/NV(I,JC)).GT.EPCON)  IDIV=IDIV+1
 42   CONTINUE
C
C------ If all densities have a small enough change, solution achieved
 47   IF(IDIV.EQ.0.AND.JTER.GT.17) WRITE(6,113) JTER
 113  FORMAT(/1X,'ITER=',I5,'   difficulty converging in odd N soln'//)
      IF(IDIV.EQ.0)  GO TO 222
C
 220  CONTINUE
C**********************************************************************

      IF(ISOL.EQ.0.AND.DCR.EQ.1) GO TO 222

      !--- Signal non-convergence here
      IF(MINJ.LE.1) WRITE(6,115) 
 115  FORMAT(' NON-CONVERGENCE in odd N solution in NORTH')
      IF(MINJ.GT.1) WRITE(6,116) 
 116  FORMAT(' NON-CONVERGENCE in odd N solution in SOUTH')
      !..WRITE(6,*) ' PT,  ALT ,   N(2D) ,    N(4S) ,    NO'
      !... print current and saved densities
      !..DO J=1,JBP1
      !..  WRITE(6,'(1X,I4,F7.1,1P,9E10.2)') J,ZR(J)
      !.. >   ,(NV(IS,J),IS=1,3),(NVS(IS,J),IS=1,3)
      !..ENDDO
      NOCON=-1    !... signals non-convergence to main
      RETURN

 222  CONTINUE
C
C------store the odd N solutions for use in vib N2 solution
      DO 88 J=1,JBP1
         N2D(J)=NV(1,J)
         N4S(J)=NV(2,J)
         IF(IPR.GE.3) NNO(J)=NV(3,J)
         IF(IPR.GE.4) O1D(J)=NV(4,J)
         IF(IPR.LT.3) NV(3,J)=NNO(J)
         IF(IPR.LT.4) NV(4,J)=O1D(J)
         JCON=J
         IF(MINJ.GT.1) JCON=MAXJ+1-J
         DO 94 I=1,IPR
           STR(I+7,JCON)=NV(I,J)
 94      CONTINUE
      IF(ZR(J).LT.0.AND.MINJ.LE.1) 
     >  WRITE(6,901) J,ZR(J),STR(8,JCON),STR(9,JCON),STR(10,JCON)
 901  FORMAT(' *CALODD* ',1X,I5,F7.1,1P,11E10.2)
 88   CONTINUE
      RETURN
      END
C::::::::::::::::::::::::  JBLK ::::::::::::::::::::::::::::::::::::::::
      BLOCK DATA JBLK
      COMMON/VAGRID/JBV,JVB,JNB
      DATA JBV, JVB, JNB /0, 0, 0/
      END
C:::::::::::::::::::::::: VERGRD ::::::::::::::::::::::::::::::::::::::::
C------- altitude grid for odd N and N2*
      SUBROUTINE VERGRD(IVLPS,GRD)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      INTEGER NION,VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      REAL GRD
      COMMON/ALT/Z(401),DT,DW,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,NION
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/VAGRID/JBV,JVB,JNB
      DATA ZUPMAX/600.0/,RE/6357./,SCL/0.0D0/,STEP/2.5/,DEN/0.0D0/
      DATA ZDOWN/90.0/  !.. Absolute lower boundary for odd N and N2v

      IF(JBV.NE.0) RETURN
      ZUP=ZUPMAX
      IF(Z(JMAX/2).LT.ZUPMAX) ZUP=Z(JMAX/2-1)
      ZR(1)=Z(1)         !.. set lower altitude boundary
      IF(ZR(1).LE.ZDOWN)  ZR(1)=ZDOWN    !.. lower boundary odd N and N2v
      !... set user specified grid steps
      IF(GRD.GE.1.0) THEN
        STEP=GRD
        IF(STEP.GE.5.0) STEP=5.0    !.. constant steps
      ENDIF      
      DH=STEP*1.0E+5       !.. alt step for equations
      !.. scales points on grid. ZUPMAX introduced 2020-07-23
      SCL=DLOG(30./STEP)/(ZUPMAX/ZR(1)-1.0)

      WRITE(17,91)   !.. note that unit 17 keeps being overwritten by FIELD
 91   FORMAT(/1X,' Vertical grid for odd N and N2(v)'/1X
     >  ,' Pt.  ALT  DELZ    BG    O1D    N2D    NO     N4S')
      JVB=1
      JNB=1
      !... loop to set up the altitude grid
      DO 930 J=1,200
         IF(ZR(J).LE.ZDOWN) JVB=J        !.. lower boundary for N2(v)
         IF(ZR(J).LE.ZDOWN) JNB=J         !.. lower boundary for odd N
         !--- grid spacing factor
         BG(J)=DEXP(-SCL*(ZR(J)/ZR(1)-1.0))
         IF(STEP.GE.5) BG(J)=1.0       !.. factor for constant steps
         DELZ=STEP/BG(J)               !.. actual altitude steps
         ZR(J+1)=ZR(J)+DELZ            !.. new altitude
         GR(J)=-980.665/(1.0+ZR(J)/RE)**2  !.. gravity

         !---- rough initial odd N profiles 
         IF(ZR(J).GE.ZDOWN) NNO(J)=2.95E+8*DEXP(-3.37E-2*ZR(J))
         IF(ZR(J).LT.ZDOWN) NNO(J)=4.54E+2*DEXP(0.1*ZR(J))
         IF(J.EQ.1.OR.ZR(J).LE.220)
     >     DEN=1.E7*EXP(-7.2E-4*(ZR(J)-200.0)**2)
         IF(ZR(J).GT.220) DEN=DEN*EXP((ZR(J+1)-ZR(J))*0.017*GR(J)/999.)
         N4S(J)=DEN
         O1D(J)=0.001*DEN
         N2D(J)=0.03*DEN

         !... tests to exit loop
         IF(ZR(J+1).GE.Z(JMAX/2)) GO TO 933  !.. exceed max field line alt
         IF(ZR(J).GE.ZUP) GO TO 933          !.. exceed max odd N alt
         IF(J.GT.VD-2) GO TO 933             !.. exceed array dimensions
         WRITE(17,'(I5,2F9.2,1p,9E10.2)') J,ZR(J),DELZ,BG(J),O1D(J)
     >     ,N2D(J),NNO(J),N4S(J)
 930   CONTINUE
 933  JBV=J-1
      RETURN
      END
C::::::::::::::::::::::::::::: VATMOS ::::::::::::::::::::::::::::::::::
      SUBROUTINE VATMOS(J250,MINJ,MAXJ,GRD,VBLAT,IC,ISOL,MAXIT,IVLPS,
     >   JBV)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      INTEGER NION,VD,J250
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      REAL AP,DEC,ETRAN,BLON,F107,F107A,SEC,D(19),T(2),CHI,SAT,GLON
     >   ,GLATD,GRD,RBLAT
      COMMON/FJS/N(4,401),TI(3,401),FF(20)
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/ALT/Z(401),DT,DW,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,NION
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      DATA SECSAV/0.0D0/

      RBLAT=VBLAT
      CALL GET_VGRID_BLAT(GRD,BLON,RBLAT,BLAT)

      !--  Get the neutral densities and interpolate ion densities
      DO 17 J=1,JBV+1
      CALL AMBS(J250,ZR(J),BLAT,D,T,CHI,SAT,GLON,GLATD)
      CHIV(J)=CHI
      HE(J)=D(1)
      ON(J)=D(2)
      N2N(J)=D(3)
      TN(J)=T(2)
      O2N(J)=D(4)
      HN(J)=D(7)
C,,,,, ion densities and te are interpolated from main prog
 15   NR(1,J)=FITDEN(1,MINJ,MAXJ,Z,N,ZR(J))
      NR(2,J)=FITDEN(2,MINJ,MAXJ,Z,N,ZR(J))
      TE(J)=FITTEM(3,MINJ,MAXJ,Z,TI,ZR(J))
      TIZ(J)=FITTEM(1,MINJ,MAXJ,Z,TI,ZR(J))
 17    CONTINUE
C
CGAG---------------------------------------------------------------|
CGAG  G. GERMANY 10/28/91 Interpolating EUV and e* production rates
      CALL VNTERP(JBV,MINJ,MAXJ)
CGAG---------------------------------------------------------------|
C---  set up production rates for odd nitrogen
      DO 35 J=1,JBV+1
         CALL PRODIT(J,1,MINJ,MAXJ)
 35   CONTINUE
      RETURN
      END
C::::::::::::::::::::::::::: CALDIF :::::::::::::::::::::::::::::::::::::
C.... This routine sets up the derivative of the flux for the solver
C.... for both odd N and N2(v). 
C.... Modified by P. Richards in April 2003 to include eddy diffusion.
      SUBROUTINE CALDIF(JBP1,IPR,AM,IC)
      IMPLICIT DOUBLE PRECISION(A-H,K-L,N-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      COMMON/DD/D(7,VD),E(7,VD),B1(7,VD),B2(7,VD),RDZ(VD),EDDY
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/ALT/Z(401),DT,DW,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,IION
      DIMENSION AM(7),EN(VD),EK(7,VD)
      DATA  BOLTZ/1.3807E-16/

      !..EDDY=0.0E6  !.. eddy diffusion coefficient (2-9E+6 cm2 s-1)

      !.. Calculate and store some things used often below
      DO 5 J=1,JBP1
        RDZ(J)=0.5*BG(J)/DH       !.. derivative factor
        EN(J)=ON(J)+N2N(J)+O2N(J) !.. total density
        DO 5 IS=1,IPR
          RMA=(ON(J)*16.+N2N(J)*28.+O2N(J)*32.)/EN(J)
          D(IS,J)=1.52E+18*SQRT(1.0/AM(IS)+1.0/RMA)*SQRT(TN(J))/EN(J)  !..&&
 5    CONTINUE

      !.. outer loop for altitude
      DO 10 J=2,JBP1
        DLTDZ=(TN(J+1)-TN(J-1))*RDZ(J)/TN(J)
        DBGDX=0.5*(BG(J+1)-BG(J-1))/DH
        D_EN_DZ=(EN(J+1)-EN(J-1))*RDZ(J)
        D2_EN_DZ2=(BG(J)/DH)**2.*(EN(J+1)-2.0*EN(J)+EN(J-1))

        !.. inner loop for species
        DO 10 IS=1,IPR
          HN2=-BOLTZ*TN(J)/AM(IS)/1.66E-24/GR(J)
          E(IS,J)=DLTDZ+1.0/HN2
          DEDZ=(BG(J)/DH)**2.*(TN(J+1)-2.0*TN(J)+TN(J-1))/TN(J)+DLTDZ*
     >      DBGDX-DLTDZ**2+((GR(J+1)-GR(J-1))*RDZ(J)/GR(J)-DLTDZ)/HN2
          DDDZ=(D(IS,J+1)-D(IS,J-1))*RDZ(J)

      !.. Eddy diffusion
          HATM=-BOLTZ*TN(J)/28.0/1.66E-24/GR(J)
          EK(IS,J)=DLTDZ+1.0/HATM
          DEKDZ=(BG(J)/DH)**2.*(TN(J+1)-2.0*TN(J)+TN(J-1))/TN(J)+DLTDZ*
     >      DBGDX-DLTDZ**2+((GR(J+1)-GR(J-1))*RDZ(J)/GR(J)-DLTDZ)/HATM
          !...DKDZ=(D(IS,J+1)-D(IS,J-1))*RDZ(J)  !.. no eddy height variation

          !.. Alternative Eddy diffusion term = d(K*N*d(n/N)/dz)/dz from
          !.. Schunk and Nagy text, where N = [N2]+[O2]+[O]
          !.. D_ED_DZ= 0   !.. derivative of eddy diffusion if alt. dependent 
          !.. ED1=D_ED_DZ - EDDY*D_EN_DZ/EN(J)
          !.. ED2=EDDY*(D_EN_DZ)**2/EN(J)**3 - 
          !.. >       (EDDY*D2_EN_DZ2+D_ED_DZ*D_EN_DZ)/EN(J)
          ED1=0  !..&&
          ED2=0  !..&&
          B1(IS,J)=DDDZ + D(IS,J) * (E(IS,J)+DBGDX) + ED1
     >       + EDDY * (EK(IS,J)+DBGDX)
          B2(IS,J)=(D(IS,J)*DEDZ+E(IS,J)*DDDZ) + ED2
     >       + EDDY*DEKDZ

          !.. if((j/10)*10.ne.j.and.j.ne.2) go to 10
          !.. write(3,99) j,zr(j),d(is,j),e(is,j),dedz,dltdz,hn2,dbgdx
 10   CONTINUE
 99   FORMAT(2X,I5,F7.0,1P,22E9.1)
      RETURN
      END
C:::::::::::::::::::::::: OFIJ :::::::::::::::::::::::::::::::::::::::
C------ this subr sets up error fns. of the t-d diff eqtns for n(2d),n(4s)
C------ no, & o(1d) densities, for  soln by the newton iterative procedure
C------ written by p. richards sept 1979
C.... Modified by P. Richards in April 2003 to include eddy diffusion.
      SUBROUTINE OFIJ(J,ILJ,IPR,JBV,IC)
      IMPLICIT DOUBLE PRECISION(A-H,K-L,N-Z)
      INTEGER VD,FEXISTS
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION PROD(6),LOSS(6),VCON(VD),FYUBV(6)
      COMMON/DD/D(7,VD),E(7,VD),B1(7,VD),B2(7,VD),RDZ(VD),EDDY
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      COMMON/ALT/Z(401),DT,DW,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,IION
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,JTER
      DATA VCON/VD*0.0/

      !-- Read in the upper boundary flux from a file
      DATA FYUBV,IR69/6*0.0,0/
      CALL TFILE(69,FEXISTS) !... see if file exists
      IF(FEXISTS.EQ.1.AND.IR69.EQ.1) THEN
        READ(69,*,ERR=177,END=177) (FYUBV(IL),IL=1,6)
        WRITE(6,*) ' FY= ',FYUBV(4)
 177    IR69 = 0
      ENDIF
 
      STS=1.0  !.. switch for steady state (may be obsolete)
      IF(ITF.GE.3) STS=0.0

      CALL PRLN(1,ILJ,J,0,IPR,PROD,LOSS,VCON,1) !.. chemistry

      !... Set up the equation for the solver
      DO 20 IS=1,IPR
        !.. upper boundary condition......*
        IF(J.EQ.JBV) NV(IS,J+1)=-E(IS,J)*NV(IS,J)/RDZ(J)
     >     +NV(IS,J-1)-FYUBV(IS)/RDZ(J)/D(IS,J)
        !.. d2n/dz2, dn/dt, and n coefficients of the differential equation
        !.. including eddy diff.
        A1=(D(IS,J)+EDDY)*
     >    (BG(J)/DH)**2.*(NV(IS,J+1)-2.*NV(IS,J)+NV(IS,J-1))
        A2=B1(IS,J)*(NV(IS,J+1)-NV(IS,J-1))*RDZ(J)
        A3=B2(IS,J)*NV(IS,J)

        !.. form error function from continuity eqtn
        F(IS)=A1+A2+A3-(NV(IS,J)-NVS(IS,J))*STS/VDT+PROD(IS)-
     >   LOSS(IS)*NV(IS,J)

        !-- write out cpts of cty eqn. the 2nd write is used after convergence
        IF((J/10)*10.NE.J.AND.J.NE.2) GO TO 20
        IF(ILJ.NE.-1.OR.JTER.NE.1) GO TO 20
        !-- a4,a5 are only used to get flux
        A4=(NV(IS,J+1)-NV(IS,J-1))*RDZ(J)
        A5=E(IS,J)*NV(IS,J)
        FLUX=-D(IS,J)*(A4+A5)
        WRITE(3,99) J,ZR(J),A1,A2,A3,A4,A5,FLUX,PROD(IS),LOSS(IS),F(IS)
 20   CONTINUE
 99   FORMAT(I4,F6.0,1P,19E11.2)
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE OMATRX(S,RHS,NEQ,IPR,MIT,JBL,JBV,IC)
C---- this program sets up the jacobian matrix for the band solver
      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
      INTEGER NEQ,VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DIMENSION RHS(1806),S(NEQ,23)
      COMMON/NNVV/F(20),NV(7,VD),NVS(7,VD)
      IBW=4*IPR-1
C
      DO 180 JZS=1,NEQ
      DO 190 KZS=1,IBW
      S(JZS,KZS)=0.0
  190 CONTINUE
  180 CONTINUE
C
      DO 80 JF=2,MIT
      J1=MAX0(2,JF-1)
      J2=MIN0(JF+1,MIT)
      DO 90 JV=J1,J2
      DO 130 IV=1,IPR
C
      L=IPR*(JF-2)
      IF(JF.LE.3)M=IPR*(JV-2)+IV
      IF(JF.GT.3)M=IPR*(JV-JF)+IV+2*IPR
      JVC=JBL+JV-1
      JFC=JBL+JF-1
C
      H=1.0E-2*NV(IV,JVC)
      NV(IV,JVC)=NV(IV,JVC)+H
      CALL OFIJ(JFC,1,IPR,JBV,IC)
      NV(IV,JVC)=NV(IV,JVC)-H
      DO 75 IS=1,IPR
      IF(JF.LE.3) S(L+IS,M)=(F(IS)-RHS(L+IS))/H
      IF(JF.GT.3) S(L+IS,M-IS)=(F(IS)-RHS(L+IS))/H
 75    CONTINUE
C
  130 CONTINUE
   90 CONTINUE
   80 CONTINUE
      RETURN
      END
C::::::::::::::::::::::::: GET_VGRID_BLAT ::::::::::::::::::::
C...... This routine returns the magnetic latitude for the vertical grid
      SUBROUTINE  GET_VGRID_BLAT(GRD,   !.. geographic latitude (in)
     >                          BLON,   !.. magnetic longitude (in)
     >                      BFREGION,   !.. magnetic latitude in F region (in)
     >                         BLATV)   !.. magnetic latitude for long GRD (out) 
      IMPLICIT NONE
      INTEGER JTIMES
      REAL GRD,BLON,GLATV,BLONV,GLAT,GLON,BFREGION,BLAT
      DOUBLE PRECISION BLATV
      DATA JTIMES/0/

       BLATV=BFREGION  !.. default location of vertical grid

      !.. see if GRD specifies a new latitude
      IF(GRD.LT.-1.00001.AND.GRD.GE.-90.0) THEN

        !.. This call finds the geographic coordinates for default lat.
        CALL BTOG(1,BLON,BFREGION*57.29578,GLAT,GLON)

        !.. establish geographic latitude for south or north
        GLATV=GRD
        IF(BFREGION.GT.0) GLATV=-GRD

        !.. get new B latitude for geographic lat. GRD
        CALL GTOB(1,BLONV,BLAT,GLATV,GLON)
        BLAT=BLAT/57.29578

      !.. Make sure new B lat is equatorward of previous lat
        IF(ABS(BLAT).LT.ABS(BFREGION)) THEN
          BLATV=BLAT  
          IF(JTIMES.LT.3) WRITE(6,90) GLATV,GLON
        ENDIF
      ENDIF

      !..WRITE(6,'(8F9.2)') GRD,BLON,BLONV,BFREGION,BLATV,GLAT,GLATV,GLON

      RETURN
 90   FORMAT(' **Vertical grid changed to (',F5.1','F5.1') geographic')
      END